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Three solution algorithms, explicit under-relaxation, point implicit, and lower-upper symmetric Gauss-Seidel, 
are Led tocomputenonequi.ibrium flow around the Apollo 4 return capsule a. the 62-km altitude pomt m Usd. 
scent trajectory By varying the Mach number, the efficiency and robustness of the solution a'gonthms were tote 
for diSLt levels of chemical stiffness. The performance of the solution algorithms degraded as the Mach number 

Gauss-Seidel algorithm deteriorates to the point that it is out performed by the point implicit method. Th 
of the viscous terms are investigated. Grid dependency questions are explored. 


Nomenclature 

= speed of sound 
= species mass fraction 
= freestream speed of sound 
= total energy 

= species heat of formation, J/kg 
= grid scaling reference length 
— molar mass of species s, gm/mole 
= heat conduction 

= universal gas constant, 8.3144 J/mole-K 
= freestream Reynolds number, Poo^oo^ref/ Poo 
= time, s 

= Cartesian velocities 
= Cartesian diffusion velocities 
= chemical source term 
= Cartesian coordinates 
= thermal conductivity 
= viscosity 
= freestream viscosity 
= density of species s 
= freestream density 
, x zz = shear stresses 

Introduction 

F UTURE human space exploration will require sending space 
vehicles to and from the moon, Mars, and beyond. The vehi- 
cles used to facilitate the return missions, whether aerobrakes or 
re-entry capsules, will be large diameter, blunt spacecraft that will 
re-enter the Earth’s atmosphere at high velocity. Tauber et al. per- 
formed trajectory studies on a 5-m-diam Mars return aerobrake that 
would enter the Earth’s atmosphere at between 12 and 16 km/s and 
experience peak heating at between 64- and 68-km altitude. Lunar 
return vehicles will experience velocities in excess of 10 km/s. These 
spacecraft will often be traveling at an angle of attack creating a true 
three-dimensional flow environment. 
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Numerical algorithms will be used to calculate the aero- and ther- 
mal loads these vehicles will encounter. The combination of large 
body diameter and high-entry velocity creates near-equilibnum flow 
conditions in some regions of the computed shock layer. These con- 
ditions lead to numerical stiffness in the governing equations affect- 
ing the stability and robustness of the numerical algorithm used to 
calculate the flowfield. The computation of three-dimensional re- 
acting flows is CPU intensive due to the complex physical models 
used and the large number of governing equations to be solved. The 
numerical algorithm must, therefore, also be efficient. 

A considerable amount of research activity has been devoted 
to developing nonequilibrium flow codes that are both robust and 
efficient. 2-6 These codes have generally been tested and applied to 
flows above 70 km or to the reproduction of ground-based experi- 
ments. . . .. f 

This study tests the stability and convergence characteristics ot 

three numerical algorithms when applied to very stiff flow condi- 
tions, a large diameter vehicle traveling at high velocity at an altitude 
of 62 km. The performance of each method is evaluated over a range 

of test conditions. 


Governing Equations 

The three-dimensional Navier-Stokes equations represent the 
conservation of mass, momentum, and energy. For chemical non- 
equilibrium flows, they include species continuity equations. The 
equations expressed in Cartesian coordinates are 
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Table 1 Specific heat curve fit constants 


Species 


a 2 

a 3 

a 4 

“5 


N 

2.4690 

1.8827c-4 

— 1 .7784f-7 

5 9797^-1 1 

-7.9413f-15 

4 9765c- 19 

O 

2.7919 

—3. 1 1 13^-4 

LI 147c-7 

— 1.5967c- 1 1 

] .1444^-15 

-4 1633c-20 

NO 

3.2943 

1 . 1040c-3 

— 3.8424f-7 

6.5532c- 11 

—5.5737c- 15 

2 425 lc-I9 

n 2 

3.3141 

1.0055c-3 

-3.7609f-7 

8.171 0^ - 1 1 

— 1.0882c- 1 4 

8 7 180c- 1 9 

O 2 

3.2489 

1.372c-3 

—5.283 le-7 

1 . 1474c-10 

— 1 .3237c- 14 

7 9561c- 1 9 

NO + 

V1 4- 

3.3277 

9. 2224c -4 

— 2.9325c-7 

4,9489^-11 

—5.0508c - 15 

3 6630c- 19 

N 2 

3.4070 

4.503c-4 

2.0053c-7 

— 8.8528c-l 1 

1.370c- 14 

1 0469c- 18 

°2 

3.3036 

1.0437e-3 

— 2.7916c-7 

8.335c- 12 

7.1209c- 15 

— 1 .01 5 1c- 1 8 

o+ 

2.5018 

— 1 .2835c-5 

2,279c-8 

— 1.391 le-l 1 

3.5489c- 1 5 

-3.7052c- 19 

N + 

2.6891 

—3. 0729c -4 

1.5156c-7 

— 3. 1635c-I 1 

3.6498c- 15 

— 2.3963c-19 

e~ 

2.5 

0.0 

0.0 

0.0 

0.0 

0.0 


a 7 


— 1.4627c -23 
6.4779e-25 
— 5.0930c-24 
-3.6331e-23 
-2.3639c-23 
-I 5729c-23 
3.9647c-23 
5. 1945c-23 
I.7067c-23 
8.321 Ic-24 
0.0 




1 .61 16c-28 
0.0 

3.930 lc-29 
5.9I53c-28 
2.7448c-28 
2.7385c-28 
-5.9278c-28 
-9.2875c-28 
— 2.8987c-28 
— 1 . 1783c-28 
0.0 


P](v + i>i)~ 


Pi(w + Wi)~ 



p n (w + w„) 

puv 

G = 

puw 

pv 2 + p 


pvw 

pvw 


pw 2 -f- p 

(e 4- p)v _ 


_ (e -J- p)w 


R — I 0, . . . , 0, x xx , r xy , x xz , q x + ujT x j — ^ p s u s h x 


S ~ ^ 0 , X xy , Xyy , Z yz , U j X yj ~ ^ p x v x h x 


T = |° 0 . r *z> hz> r zz< qz + UjX zj - PsW x k x 

W = [w { U) ny 0,0, 0,0] 

In the preceding expressions 

UjXxj — UT XX + vx xy + w X xz 

Expressions for the shear stresses are provided in the Appendix. 

The chemical source terms w x represent the production of species 
from finite rate chemical reactions. In this study, a seven-species air 
chemistry model is used. This was chosen over an 1 1 -species model 
in an effort to reduce computer time requirements. Charge neutrality 
is assumed throughout the flow, so only six species, (N, O, NO, N 2 , 
0 2 , NO + ), are included in the equation set. Six chemical reactions 
are evaluated to determine the species production terms. They are 

N 2 -hM ^N+N-f M 

o 2 + a/^o + o+ a/ 

NOf O + M 

NO + O ^ N 0 2 
N 2 + O ^ N + NO 
N-fO^ NO+ + e~ 

The quantity M in the expressions can be any one of the species. 
Forward and backward reaction rates for the reaction set are obtained 
from Park. 7 

Pressure is obtained from the expression 



The equation relating temperature to total energy is 

e ~ ^2 R j dT + Psh" + \p{u 2 + v 2 + w 2 ) (2) 

Values of species specific heats are obtained by curve fit re- 
lations. Previously published specific heat curve fits 8 ' 9 use a seg- 
mented approach; the curve fit is divided into temperature bands, 
each band employing a different mathematical relation to obtain the 
specific heat. Segmented curve fits can have discontinuities at the 
segment boundaries, and conditional relations must be incorporated 
into the code to identify which mathematical relation should be used. 

An effort was undertaken to develop new curve fit relations that 
would provide accurate values of species specific heat with no seg- 
mentation of the curve fits. One mathematical expression is used 
throughout the entire temperature range, from 100 to 20,000 K. Un- 
der chemical nonequilibrium conditions, postshock temperatures 
will not exceed 20,000 K for any realistic Earth entry velocity. The 
curve relations take the form 

C r /R = a { +a 2 T -ha^T 2 +a 4 T y +a 5 T 4 + +a 7 r 6 -\-a^T 1 (3) 

The constants were obtained by matching reference values for 
C p / R at eight temperatures from 0 to 20,000 K. Below 6000 K the 
reference values were taken from JANAF thermochemical data. 10 
Reference values for the neutral diatomic species, N 2 , 0 2 , and NO, 
above 6000 K were calculated by Jaffe, 11 using rigorous quantum 
mechanical definitions of the molecular partition functions. Ref- 
erence values tor the atomic and ionized diatomic species above 
6000 K were taken from Ref. 8. 

Because no segmentation is used, the specific heat curve is con- 
tinuous throughout the entire temperature range. Species enthalpy 
or internal energy can be obtained by integrating Eq. (3). Values of 
the constants, a\ , a 2 for 1 1 -species air are listed in Table 1. 

Nondimensional Numbers 

Two nondimensional parameters are used in the analyses pre- 
sented in this paper. The first is the Damkohler number defined as 

D (l ~ T lraa(i /r chcm (4) 

The quantity T trans is the characteristic time of translation. It can 
be evaluated for each computational cell by dividing the length 
of the cell by the mean flow velocity. The parameter r C h em is the 
characteristic time of chemical relaxation. This can be computed 
for each chemical reaction and is a function of total density, species 
density, and temperature. Two limiting cases of Damkohler number 
are if D a 1 , the flow can be considered frozen, and if D u 1, 
the flow is in chemical equilibrium. 

The second nondimensional number used in this study is the 
Courant— Friedrichs-Lewy (CFL) number defined as 

CFL = |A ma JAt/A.r 

where A max is the maximum of the eigenvalues. Explicit methods 
cannot run at a time step that yields a CFL number greater than one, 
Implicit methods can run at CFL numbers greater than one. 
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Solution Algorithms 

The governing equations are transformed from Cartesian (x , y, z) 
to a generalized 17, C) coordinate system. The generalized 
Navier-Stokes equations are 

i)Q DE i)F 3G _ J_ 

l)t + 3£ + <h\ + 3C Re 

The transformed Q and W vectors a 

Q = Q/J 


HR as ar 

— “i - — 

at a n a< 

re given by 
W = W/J 


+ w 


where 

r x . . aw' 

L = / + A/ [r fl + 0, + r t .]/ — A-_ j — — C t-1 — — 

£> = [/ + At ([r„ + r,, + r, .]/)]-’ 

U = I + At[[r„ + r h + r,-]I + A~ +1 + fl~ + , + C t 7 +1 ] 

The approximate split flux Jacobians A + and A are computed 
from the equation 


and the transformed flux vectors are 


^■ E+ 5 -r f ‘ c 

n S + ZzT 

R ~ J 

ri.E + t] v F + r) z G 
F ~ J 

„ r) x R + ri y S+r) z f 
S ~ J 

G ~ J 

ft it + SyS + SlT 
1 - J 


The quantity J in the expressions is the Jacobian of the coordinate 
transformation. 

The simplest solution algorithm is the explicit Euler algorithm. 
This method evaluates the flux and source term vectors at the cur- 
rent, or n, time step. The three-dimensional explicit Euler solution 
algorithm written in generalized coordinates is 

SQ = -At[S !i E+& n F+8 < G-W~(\/Re)(8 i :R + 8 n S + 8 l; T)] (5) 

where 

&Q = Q" + ' - O' 


A* = (A±r„/)/2 (8) 

where r a is the spectral radius equal to the absolute value ot the 
largest eigenvalue, 

r a = \U\+c^f+^f+^ (9) 

In Eq. (10), U is the contravariant velocity in the £ direction, 
% x u + % y v + £ z ui. Similar expressions are used to calculate the B± 
and C± matrices. To obtain 8Q, the [L] and [U] elements are solved 
by sweeping from one comer of the computational domain to the 
other. One matrix inversion is required per point per step because 
the (3W/3g) matrix must be inverted. 

The third method tested was developed by Palmer. 14 This tech- 
nique maintains the stability of a nonequilibrium algorithm within 
the context of the explicit formulation . Thi s explicit under-relaxation 
algorithm evaluates the species density equations as well as the total 
density conservation equation explicitly. From these equations the 
changes in species mass fractions are calculated 

&Ci = (8pi-Ci&p)/p 00 ) 


The term 8^ indicates a spatial difference in the £ direction. If no 
source term is present, this method has a maximum allowable CFL 
number, based on the fluid dynamics, of one. 

For nonequilibrium flows, the species production terms obtained 
from finite rate chemical reactions introduce an additional stiffness 
into the equation set causing the algorithm to be unstable except for 
very small CFL numbers. For this reason, the explicit Euler solution 
algorithm is widely regarded as unsuitable for the computation of 
nonequilibrium flows. 

To overcome the time step restrictions imposed by the chemical 
source terms, the source terms can be evaluated implicitly, at the 
n + 1 time step. A Taylor series expansion is performed on this term 
and the resulting algorithm, called point or chemistry implicit, is 
given by 


/ 


aw 

ZQ 


8Q = RHS 


( 6 ) 


where the right-hand side (RHS) in Eq. (6) contains the explicit 
terms, the elements on the right-hand side of Eq. (5). In this study 
the inviscid flux differences are evaluated using Van Leer flux vec- 
tor splitting. 12 Full viscous terms are included as well as a binary 
diffusion model. 

This algorithm is called the point implicit method because the 
implicit terms are not spatially differenced but are evaluated point 
by point. Because the chemistry terms are handled implicitly, this 
algorithm has a CFL number limitation of one. 

To achieve CFL numbers greater than one, it is necessary to 
evaluate the fluid dynamic flux terms implicitly. One method that 
has gained popularity in recent years is the lower-upper sym- 
metric Gauss-Seidel (LUSGS). It was developed by Yoon and 
Jameson. 13 The method splits the Jacobian matrices, A = (OE/flg), 
B = (.'>F/ag), and C = (3G/3g), into positive and negative com- 
ponents and factors the resulting split matrices into three submatri- 
ces. The resulting algorithm is given by 


(7) 


If the absolute value of the maximum change in species mass fraction 
is greater than a prescribed value tol then the changes in species mass 
fraction are scaled 


8c i = 


8c x 

|<5c max I 


* tol 


(ID 


The parameter tol is an under-relaxation factor with a value between 
0 and 1 . Additionally, no species mass fraction is allowed to become 
negative. Species densities are updated by 

=c? + V +1 02) 


In effect, this technique reduces the chemistry time step in regions 
of stiffness and has been applied to a variety of hypersonic flow 
computations. 6, 14,15 


Boundary Conditions 

The boundary conditions used in the calculations to come were 
as follows: Along the inflow (k = k max) plane, freestream values 
are maintained. Symmetry conditions are used for the j = 1 and 
j - j max planes. Along the outflow (i = * max) plane and along 
the singular line (/ = 1), values arc obtained by extrapolation. A 
constant temperature of 1 500 K was maintained on the body surface 
that was assumed to be noncatalytic. Nonslip and zero pressure 
gradient conditions were enforced. 

Implicit boundary conditions are incorporated into the LUSGS 
option. At the body surface, they are based on the relations 

&P\ = &p 2 8T[ — 0 8c x \ -- 8c s 2 

8pu\ = 8pv i = &pw\ = 0 

The subscript 1 in the preceding relations refers to the body surtace 
and the 2 to the next interior plane. The expressions are used to 
obtain values of 8Q at the body surface. 


LDU8Q = RHS 
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Results 

Specific Heat Data 

Results generated using the new curve fit relations are com- 
pared against two previously published segmented curve fits. 
Balakrishnan* developed curve fit relations up to 50,000 K using 
from five to seven temperature bands. The curve fits have the form 

C p — a + bT + c/T 2 ( 13 ) 

The second curve fit relations, used by Gnofto et al., y have the form 

Cp/R — a \ + QiT + a^T 2 -f- a 4 T* + a 5 T 4 (14) 

and are used to calculate species specific heats up to 35,000 K using 
five temperature bands. 

Figure la compares values of the specific heat of atomic nitro- 
gen calculated using the three curve fit relations. The segmented 
curve fits exhibit discontinuities at the temperature band interfaces, 
particularly at 6,000 and 10,000 K. The nonsegmented curve fit is 
continuous throughout the entire temperature range and reproduces 
the JANAF specific heat data with a maximum error of 2.2%. At 
temperatures above 6,000 K, the three curve fits give similar data, 
except for the Balakrishnan curve fit in the temperature range of 
10,000-15,000 K. 

Figure lb shows specific heat data for NO. Again the segmented 
curve fits show discontinuities at the 3,000, 6,000, and 10,000 K 
temperature band segment interfaces. Above 12,000 K only the un- 
segmented curve fit relations match the reference values of Jaffe, 11 
which include effects such as a nonrigid rotor rotational model and 
vibration-rotation coupling. Vibration-rotation coupling produces 
a significant negative contribution to the specific heat of nitric oxide 
above 12,000 K. The nonsegmented curve fit values are in good 
agreement with the reference values along the entire temperature 
range. 




Fig. 1 Specific heat comparison: a) atomic nitrogen and b) nitric oxide. 


Solution Algorithm Comparisons 

The solution algorithms were tested by computing flow over the 
base of the Apollo 4 return capsule at the 62-km-altitude point 
along its descent trajectory. The capsule geometry is the intersec- 
tion of a 4.69-m-radius sphere with a 33-deg half-angle cone. 15 The 
cone-sphere intersection is rounded into a 19.56-cm-radius shoul- 
der. The freestream conditions at 62 km are given as 

r„ = 4.69 m 

Po c = 2.407 x 10 -4 kg/m 3 

Poo = 16.69-^ 
m 2 

T x = 241.5 K 
c x = 311.5 m/s 
mass fraction N 2 = 0.7656 
mass fraction 0 2 = 0.2344 

By varying the freestream Mach number, it is possible to simulate 
different levels of chemical stiffness. Table 2 displays the computed 
postshock stagnation line Damkohler number at 62 km for each of 
the six chemical reactions at Mach numbers of 15, 30, and 40. Mach 
40 corresponds to a freestream velocity of 1 2.5 km/s. For most of the 
reactions, the Damkohler number increases with increasing Mach 
number. Increasing Damkohler number indicates increasing chem- 
ical stiffness. The ionization of nitric oxide reaction is the stiffest 
reaction at Mach 30. The Damkohler number of the thermal dissocia- 
tion of nitric oxide reaction, reaction 3, increases more than an order 
of magnitude with each increased Mach number and is the stiffest 
reaction at Mach 40. There is a large increase in Damkohler number 
for reactions 1,2,3, and 6 when the Mach number is increased from 
30 to 40. This is in part due to the fact that the seven-species model is 
inadequate at 12.5 km/s and leads to unrealistically high postshock 
temperatures. 

The 49 x 11 x 49 grid used in the first set of parametric studies 
is shown in Fig. 2. Only every third grid line in the body normal 
and streamwise directions are shown for figure clarity. Calculations 


Table 2 Postshock stagnation line Damkohler number, 62 km 



Mach = 1 5 

Mach = 30 

Mach = 40 

n 2 +m=n+n+m 

6.8 x 10~ 4 

0.38 

97.9 

o 2 +m=o+o+m 

5.1 

67.2 

363 

N0+M=N+O+M 

0.82 

246 

10,180 

N0+0=N+C>2 

4.1 

103 

44.7 

n 2 +o=n+no 

43.7 

67 

18.4 

N+0=NO + +e ~ 

911 

463 

989 



Fig. 2 49 x 11 x 49 grid. 
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Fig. 3 Convergence histories, point implicit method: a) Mach 15, 
b) Mach 30, and c) Mach 40. 

performed using this grid were performed at zero angle of attack so 
only a 90-deg section of the Apollo vehicle is used. 

Figure 3 presents convergence histories of the point implicit 
method at Mach 15, 30, and 40. In these and subsequent figures, 
the maximum level of residual for each calculation was scaled to 
unity to aid the comparison at different Mach numbers. Figure 3a 
shows the results at Mach 15. The point implicit solution converges 
in a smooth manner while running at maximum CFL numbers of 0.8 
and 0.9. The L2 norm of the energy residual, defined as the square 
root of the sum of the energy residual at each point in the compu- 
tational domain squared, converges eight orders of magnitude after 
running 16,000 steps at CFL = 0.9. 


At Mach 30, the stiffness of the chemistry terms increases. As 
shown in Fig. 3b, the convergence rate at a CFL number of 0.9 
slows after a three order of magnitude drop in residual. The rate 
of convergence then increases but slows again after converging an 
additional four orders of magnitude. Lowering the CFL number did 
not eliminate the decrease in the rate of convergence. This same 
behavior is apparent when the maximum CFL number was reduced 
to 0.8 and 0.59. The overall convergence rate of the point implicit 
method is slightly less than that seen at Mach 15. 

Results when the point implicit method was run at Mach 40 are 
shown in Fig. 3c. The performance of the method is similar to that 
seen at Mach number 30 except the second decrease in the rate of 
convergence occurs at a later point in the computation. 

The performance of the explicit under-relaxation method is ex- 
amined in Fig. 4. The value of the under-relaxation parameter tol 
seen in Eq. (1 1) was set to 0.01 . The convergence histories for Mach 
1 5 are shown in Fig. 4a. The solution converges in a smooth manner 
at maximum CFL numbers of 0.55, 0.68, and 0.8, but when the CFL 
number is raised to 0.9 the convergence profile flattens out after 
dropping just over three orders of magnitude. 

When the Mach number is raised to 30 the convergence difficulties 
occur even when the CFL number is reduced to 0.27. This is shown 
in Fig. 4b. The oscillating behavior of the residual is clearly evident 
at CFL numbers of 0.6 and 0.8. The same trend was seen at Mach 
40. For no value of CFL number would the solution converge. 

The explicit under-relaxation algorithm is an approximate 
method. The under-relaxation parameter tol has the same effect 
as the implicit chemistry matrix does in the point implicit algo- 
rithm. Evaluating the species density updates using tol is not as 
rigorous or as time consuming as filling and inverting the implicit 
chemistry matrix. When the chemistry is not very stiff this approx- 
imate nature does not adversely affect the solution process. When 
the chemistry becomes stiff the explicit under-relaxation method 
breaks down. What happens is that the code over- or undershoots 
the proper chemical state of the gas, and the solution oscillates be- 
tween two chemical states. This happens throughout the flowfield 
but particularly in regions of high chemical stiffness. The explicit 
under-relaxation method never converges to a single steady-state 
solution. 

A parametric study on the under-relaxation parameter tol was 
performed at Mach 30 to see if the value of the parameter affected 
the convergence rate of the method. Figure 4c shows the conver- 
gence histories obtained using three values of the under-relaxation 
parameter. All three residuals flatten out and begin an oscillatory 
behavior. The solution using the lowest value of tol, 0.00 1, con- 
verges nearly one order of magnitude further before flattening out, 
but none of the tol values permitted full convergence to occur. 

Comparing Figs. 3 and 4, both the point implicit and the explicit 
under-relaxation methods experience a leveling off of the rate of 
convergence at Mach 30. The point implicit is able to overcome this 
after some time and then is able to continue to converge. The explicit 
under- relaxation cannot, and the solution from this point oscillates. 

The performance of the LUSGS algorithm is presented in Fig. 5. 
Figure 5a shows the convergence histories at Mach 15. At maximum 
CFL numbers of 20, 40, and 90, smooth convergence is achieved. 
The method became unstable, however, when the CFL number was 
increased to 180. 

Figure 5b shows the convergence histories at Mach 30. The resid- 
ual curve now flattens after converging three orders of magnitude. 
The residual then drops another three orders of magnitude at this 
reduced rate before dropping rapidly for about 300 iterations after 
which the residual curve flattens out again. The performance of the 
LUSGS method diminished when the Mach number was increased 
from 15 to 30. At Mach 30, 50% more iterations were required to 
achieve the level of convergence obtained at Mach 15. Lowering the 
CFL number did not affect the shape of the residual profile. 

The performance of the LUSGS method degrades severely at 
Mach 40, shown in Fig. 5c. When the algorithm was run at a CFL 
number of 23, the residual profile leveled off after converging two 
orders of magnitude, similar to the behavior seen with the explicit 
under-relaxation method. The convergence problem remained when 
the CFL number was reduced to 5.8. Only when the CFL number 
was reduced to 2.3 did the performance improve. At a CFL number 
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Fig. 4 Convergence histories, explicit under-relaxation method: 
a) Mach 15, b) Mach 30, and c) constant parametric study. 



Fig. 5 Convergence histories, LUSGS method: a) Mach 15, b) Mach 
30, and c) Mach 40. 


ot 2.3, the residual profile was similar to that produced by the point 
implicit method. 

There are several possible reasons besides chemical stiffness why 
the solution algorithm performance worsened w hen the Mach num- 
ber was increased. In the LUSGS method, the implicit flux Jacobians 
of the inviscid fluxes are split into positive and negative compo- 
nents and then differenced using one-sided differences. The viscous, 
which are centrally differenced, are not treated implicitly. The vis- 
cous terms might, therefore, affect the performance of the method. 

A comparison was made in the behavior of the LUSGS method 
with and without viscous terms. The algorithm was run at Mach 30 
at CFL numbers of 36 and 90. The results are shown in Fig, 6a. 


There is very little difference in the residual profiles with or without 
viscous terms, so viscous terms would not appear to be the cause 
of the diminished performance. They do have some effect on the 
stability of the method. When the viscous terms were removed, it 
was possible to run the code at a CFL number of 179. The results 
obtained without viscous terms at CFL numbers 36, 90, and 179 are 
shown in Fig. 6b. There is not a substantial difference in the three 
residual histories. 

Another possible reason for the diminished performance could be 
grid resolution effects. The spatial resolution might be insufficient 
to resolve the chemical gradients. Computations were performed 
using an axisymmetric version of the code with a 49 x 49 grid. 
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Fig. 6 Effect of viscosity, LUSGS method: a) viscous effects and b) 
inviscid LUSGS. 

equivalent to one plane of the three-dimensional grid. The compu- 
tation was repeated using an 89 x 89 grid. This represents a uniform 
increase in the spatial resolution throughout the computational do- 
main. An axisymmetric computation was performed to reduce the 
computer time necessary, but the axisymmetric solution was found 
to converge in a similar fashion to a three-dimensional zero-angle- 
of-attack calculation. Doubling the number of grid points reduces 
the Damkohler number in the flowfield by reducing the characteris- 
tic time of translation in each computational cell. 

Figure 7 shows the results of this test using the point implicit and 
LUSGS methods. Solutions using the 49 x 49 and 89 x 89 grids 
show similar convergence profile shapes. The 49 x 49 and 89 x 89 
grid solutions for both methods show a leveling off in the rate of 
convergence after the residual drops three orders of magnitude. The 
point in the computation when the residual levels off with the 89 x 89 
grid occurs after twice the number of iterations as the 49 x 49 grid. 
Increasing the spatial resolution by a factor of four everywhere in the 
computational domain does not seem to alleviate the performance 
dropoff of either method. The uniformly increased spatial resolution 
reduces the convergence rate while maintaining the same profile 
shape for both methods. 

Another possibility is that grid dependencies in the shock region 
are affecting the convergence characteristics. The highest temper- 
atures and, therefore, fastest reaction rates occur behind the bow 
shock. The adaptive grid code SAGE 17 was used to increase the 
spatial resolution in the shock region by adapting the original grid 
to the solution. This reduces the Damkohler number in the region 
of the bow shock. Results using the explicit under-relaxation and 
LUSGS methods with the original and adapted grids at Mach 30 are 
shown in Fig. 8a. There is no improvement in the performance of the 
explicit under-relaxation method when an adapted grid is used. The 



Fig. 7 Grid dependency study, Mach 30. 




Fig. 8 Effect of grid adaption: a) grid comparison, Mach 30 and b) 
grid comparison, LUSGS, Mach 40. 

initial convergence rate is lower, and both the original and adapted 
grid residual profiles asymptote to the same value. The convergence 
profile for the LUSGS method using the adapted grid is smoother 
than that obtained using the original. Both solutions obtain the same 
overall rate of convergence. 

The effect of using the LUSGS method with an adapted grid at 
Mach 40 is shown in Fig. 8b. The solution does not converge using 
the adapted grid when the algorithm is run at a CFL number of 55 
or 21 . When the CFL number is reduced to 5.5, the performance of 
the method is improved using the adapted grid, although the overall 
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Fig. 9 Algorithm comparison in terms of computer time usage' a) 
Mach 15, b) Mach 30, and c) Mach 40. 

convergence level achieved after 3000 steps is not significantly dif- 
ferent from that obtained using the original grid at a CFL number 
of 2.3. 

Figure 9 compares the solution algorithms against each other 
in terms of convergence vs Cray C-90 CPU time expended. The 
results presented here might vary from code to code depending on 
the level of vectorization employed. The point implicit and LUSGS 
algorithms both required one 10x10 matrix inversion per point per 
step. The matrix inverter was highly vectorized, utilizing an inner 
loop of length 26,41 1 (49 x 1 1 x 49). All three methods used the 
same evaluation of the explicit terms, the terms on the right-hand 
side of Eq. (5). In this study, the explicit under-relaxation method 


required 72.9 /xs per point per step or 1 .44 s per step on a 49 x 11 x 49 
grid. The point implicit method used 1 .93 s per step, and the LUSGS 
algorithm required 5.07 s per step. 

Results at Mach 15 are presented in Fig. 9a. The LUSGS algo- 
rithm produces an eight order of magnitude drop in the L2 norm of 
the energy residual in one-third the CPU time of the point implicit or 
explicit under-relaxation methods. The explicit under-relaxation al- 
gorithm outperforms the point implicit method by a smaller margin. 

At Mach 30, the performance of the LUSGS method has deteri- 
orated, but the algorithm still produces an eight order of magnitude 
drop in the residual in half the CPU time of the point implicit method 
as shown in Fig. 9b. The explicit under-relaxation method failed to 
produce a converged solution at this Mach number. At Mach 40, the 
LUSGS algorithm no longer outperforms the point implicit method. 
To achieve a converged solution, it was necessary to severely reduce 
the CFL number with the LUSGS method. From Fig. 9c, it is not 
apparent whether the LUSGS method would achieve the seven or- 
der of magnitude drop in residual achieved by the point implicit 
algorithm. 

The three solution algorithms yielded essentially the same 
steady-state solution, as they should since all three use the same 
right-hand-side elements. When the explicit under-relaxation algo- 
rithm would break down, the solution would oscillate between two 
chemical states. At a given point the temperature might oscillate 
between, say, 5200 and 5220 K with a corresponding oscillation in 
species mass fraction and the other flow quantities. 


Summary and Conclusions 

Three solution algorithms are used to compute nonequilibrium 
flow around the Apollo 4 return capsule at 62-km altitude. By vary- 
ing the Mach number, the efficiency and robustness of the solution 
algorithms were tested for different levels of chemical stiffness. 

The point implicit method was able to converge at all Mach num- 
bers at CFL numbers close to the explicit limit of 1.0. The shape 
of the residual profiles was not affected by increasing CFL number. 
The method showed only a slight performance degradation with 
increasing chemical stiffness. 

The explicit under-relaxation algorithm achieved convergence 
only at Mach 1 5 and then only for a CFL number below 0.9. Above 
Mach 15, the method would only converge one order of magnitude 
before residua] profile began oscillating. 

Tlie CFL number does not affect the shape of the LUSGS residual 
profile at Mach 30 and below. Effective convergence is possible at 
high CFL numbers. At Mach 40, the LUSGS method only converged 
if the CFL number was reduced to a low value, 2.3 in this study. 
Even then the residual decreases only four orders of magnitude 
before leveling off. 

The fact that the viscous terms were evaluated explicitly did not 
affect the convergence rate or account for the performance dropoff 
of the LUSGS method at Mach 30 although it did affect the stability. 
When the solution algorithm is run inviscid it is possible to use a 
higher CFL number. 

The performance trends of the three solution algorithms seem 
to be grid independent. When the spatial resolution was uniformly 
increased by doubling the number of grid points in each direction, 
the convergence rate of the point implicit and LUSGS algorithms 
slowed, but the shape of the residual profiles remained unchanged 
The use of a solution-adapted grid did not improve the performance 
of the explicit method. The residual sill flattened out after converging 
one order of magnitude. For the LUSGS algorithm the residual pro 
file was smoother using the adapted grid, but the overall convergence 
rate was unchanged at Mach 30. At Mach 40 the LUSGS method 
showed a slight performance improvement using the adapted grid 

” ut J h T e T ^ e o ° f COnVCrgencc was stiI1 much lower than that seen with 
the LUSGS algorithm at Mach 30. 

At low levels of chemical stiffness, the LUSGS algorithm outper- 
forms both the point implicit and explicit under-relaxation meth- 
ods. At Mach IS, LUSGS used one-third the Cray C90 computer 
time to achieve the same level of convergence. As chemical stiff- 
al1 ,hree methods show a performance dropoff. At 
ac the explicit under-relaxation method cannot produce a 
converged solution. The convergence rate of the LUSGS method 
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slows more rapidly than the point implicit, but at Mach 30 LUSGS 
can still achieve a given level of convergence in half the computer 
time. 

As the Damkohler numbers continue to increase, the performance 
of the LUSGS method continues to degrade relative to point implicit. 
At Mach 40 there was no advantage to using LUSGS. The rate of 
convergence of the LUSGS algorithm was slightly less than that 
of point implicit, and LUSGS could not achieve the same level ot 
convergence. 

If an 1 1 -species air model were used, the maximum temperature 
level would decrease. This would tend to reduce chemical stiff- 
ness in these high-temperature regions. It is unknown whether the 
overall chemical stiffness would be reduced without evaluating the 
Damkohler number of the reactions involving the additional ion- 
ized species. The generation and inversion of the implicit matri- 
ces would be significantly more expensive in terms of computer 
time. 


Appendix: Viscous Terms 

The shear stresses are defined as 
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Expressions for the heat conduction terms are 
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